home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / autoseq.src / autoseq.C < prev    next >
Text File  |  1996-07-05  |  48KB  |  1,650 lines
  1.  
  2. #undef WIN_MAC
  3.  
  4. //    ============================================================================
  5. //    autoseq.C                                                        80 columns
  6. //    Reece Hart    (reece@ibc.wustl.edu)                                tab=4 spaces
  7. //    Washington University School of Medicine, St. Louis, Missouri
  8. //    This source is hereby released to the public domain.  Bug reports, code
  9. //    contributions, and suggestions are appreciated (to email address above).
  10. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  11. //    Takes a single filename of a tracefile and outputs individual files of the
  12. //    components of the tracefile.  Any format supported by CTraceFile is
  13. //    supported here.
  14. //
  15. //    autoseq was intended to be a moderately robust and complete testbed in which
  16. //    this project could be designed, implemented, and tested.  It was not
  17. //    designed to be everyone's solution to chromatogram analysis.  Nonetheless,
  18. //    it's probably satisfactory for users who seek access to specific elements
  19. //    of tracefiles and those who wish to implement their own peak picking
  20. //    routines.
  21. //
  22. //    SLATED IMPROVEMENTS
  23. //    *    AUTOSEQ env. variable for processing
  24. //
  25. //    MODIFICATION HISTORY
  26. //    93.11.11    Reece    First release
  27. //    ========================================|===================================
  28.  
  29. #include    "CTraceFile.H"
  30. #ifdef WIN_MAC
  31. #include    "CTrace.C"  // dgg
  32. #endif
  33. #include    "DNA.H"
  34. #include    "FileFormat.H"
  35. #include    "RInclude.H"
  36. #include    "RInlines.H"
  37. #include    "Definitions.H"
  38. #include    <iomanip.h>
  39. #include    <iostream.h>
  40. #include    <stdio.h>
  41. #include    <stdlib.h>
  42. #include    <string.h>
  43.  
  44. const char* VERSION        = "94.05.03";
  45. const char* COMMENT_PREFIX = "// ";            // start of line containing comment
  46. const short    FNLEN        = 200;                // max filename length
  47. const short SUFXLEN        = 5;                // max suffix length
  48. const char    BASE_TOKEN    = '#';                // placeholder for base letter in fn
  49. const char    HYPHEN        = '-';                // command line switch marker
  50. const int    DIGEST_NO_PROC = -1;            // returned if we shouldn't process
  51. const tracepos START_AT_PRIMER = -1;        // set leftCutoff to primer position
  52.  
  53.                                             // command line switches:
  54. const char* basesw        = "b";                // base selector
  55. const char* blsw        = "bl";                // baseline
  56. const char* bpdsw        = "bpd";            // base position deltas
  57. const char* bpsw        = "bp";                // base positions in trace
  58. const char* d1sw        = "d1";                // 1st deriv
  59. const char* d2sw        = "d2";                // 2nd deriv
  60. const char* fmtsw        = "fmt";            // format switch
  61. const char* fssw        = "fs";                // file sequence (stored in file)
  62. const char* ftsw        = "ft";                // fluorescence trace
  63. const char* fwsw        = "fw";                // fluorescence model breadth
  64. const char* headsw        = "h";                // include headers with output
  65. const char* helpsw        = "help";            // cl switches for
  66. const char* indsw        = "i";                // individual traces
  67. const char* leftsw        = "l";                // left cutoff
  68. const char* pmcsw        = "pmc";            // peak mean coeff for threshold
  69. const char* outprfxsw    = "o";                // output prefix
  70. const char* prunesw        = "p";                // prune
  71. const char* pssw        = "ps";                // predicted seq (pred by CPeakList)
  72. const char* ptsw        = "pt";                // peak trace
  73. const char* quietsw        = "q";                // quiet (opposite of verbose)
  74. const char* rightsw        = "r";                // right cutoff
  75. const char* rtsw        = "rt";                // residual trace
  76. const char* scalesw        = "s";                // scale switch
  77. const char* smoothsw    = "sm";                // smooth switch
  78. const char* tfssw        = "tfs";            // tracefile summary
  79. const char* tssw        = "ts";                // trace summary
  80. const char* versionsw    = "version";        // get version info
  81. const char* vsw            = "v";                // verbose
  82. const char* xformsw        = "x";                // transpose
  83. const char* zerosw        = "z";                // zero threshold
  84.  
  85. const tracepos    FLUOR_WIDTH_DEFAULT    = 50;    // exp. flrsnc Gaussian is this wide
  86. const tracepos    SEPARATION_DEFAULT    = 4;    // est'd lower limit of separation
  87. const double    PMC_DEFAULT            = 1.0;    // peak threshold = PMC_DEFAULT*avg
  88. const double    Z_THRESH_DEFAULT    = 0.0;    // 1st-deriv zero epsilon
  89. const ushort    MAX_NBHD            = 2;    // default neighborhood for deltas
  90. const short        SMOOTH_ITER_DEFAULT = 2;    // # of times to smooth
  91.  
  92. const double SCALES[NUM_BASES] =            // magic scaling values determined
  93.     {                                        // empirically
  94.     17296.90264180267,                        // A
  95.     22057.42514438165,                        // C
  96.     16231.00021563674,                        // G
  97.     7007.431127025081                        // T
  98.     };
  99.  
  100. const xform_mtx    ORTHO_MATRIX =                // default 4x4 orthogonalization
  101.                                             // matrix from David States.
  102.     {                                        // see Transform() for format
  103.     //                CONTRIBUTION TO
  104.     //         [,A]     [,C]     [,G]     [,T] 
  105.         {     20.0,     -9.0,     -4.0,      0.0    },    // FROM    [A,]
  106.         {     -7.0,     10.0,      0.5,      0.0    },    //         [C,]
  107.         {    -11.0,      0.0,     15.0,     -2.0    },    //         [G,]
  108.         {      0.0,      0.0,    -12.0,     14.0    }    //        [T,]
  109.     };
  110.  
  111. struct    opts_t;                        // structure for options passed to autoseq
  112.  
  113. //
  114. // Prototypes
  115. //
  116. void    help(void);
  117. int        process(const char *FN, opts_t& options);
  118. void    Fatal(char* error)
  119.             { cerr << "autoseq: FATAL: " << error << endl; exit(1); }
  120. void    Warning(char* error);
  121. void    Warning(char* error)
  122.             { cerr << "autoseq: WARNING: " << error << endl; }
  123. inline
  124. ostream& Banner(ostream& os)
  125.     {
  126.     CTraceFile        tf;
  127.     CTrace<int>        t(0);
  128.     CSequence<int>    s;
  129.     CPeakList        pl;
  130.     os<<COMMENT_PREFIX<<"autoseq       version "<<VERSION<<endl;
  131.     os<<COMMENT_PREFIX<<"CTraceFile   version "<<tf.Version()<<endl;
  132.     os<<COMMENT_PREFIX<<"CTrace       version "<<"unspecified"<<endl;
  133.                                              //<<t.Version()<<endl;
  134.     os<<COMMENT_PREFIX<<"CSequence    version "<<s.Version()<<endl;
  135.     os<<COMMENT_PREFIX<<"CPeakList    version "<<pl.Version()<<endl;
  136.     os<<COMMENT_PREFIX<<"by Reece Hart (reece@ibc.wustl.edu)"<<endl;
  137.     os<<COMMENT_PREFIX<<"this executable and its source are public domain.\n";
  138.     return os;
  139.     }
  140.  
  141. struct opts_t
  142.     {
  143.     // NOTE: bases uses bits 0-3 for bases A,C,G,T respectively so that
  144.     // X=ACGT from DNA.H enum is the corresponding bit and 0x01<<X is a bitmask
  145.     // for a particular base
  146.     bool    verbose;
  147.     bool    firstActionSwitch;
  148.     format_t format;
  149.     bool    headers;
  150.     char    prefix[FNLEN];
  151.     bool    pickPeaks;
  152.  
  153.     short    bases;
  154.     tracepos left;
  155.     tracepos right;
  156.  
  157.     // modifiers determine data processing before picking & how they're picked
  158.     bool    bl;            double    baseline;
  159.     bool    smooth;                                    short    smooth_iterations;
  160.     bool    scale;        char scaleFN[FNLEN];        double    scales[NUM_BASES];
  161.     bool    xform;        char xformFN[FNLEN];        xform_mtx matrix;
  162.     double    PMC;
  163.     double    zthresh;
  164.  
  165.     // operations which cause output
  166.     // at least one of the following must be true
  167.     bool    bp;            char bpFN[FNLEN];
  168.     bool    ind;        char indFN[FNLEN];
  169.     bool    d1;            char d1FN[FNLEN];
  170.     bool    d2;            char d2FN[FNLEN];
  171.     bool    fs;            char fsFN[FNLEN];
  172.     // these operations cause peaks to be picked:
  173.     bool    ts;            char tsFN[FNLEN];
  174.     bool    pt;            char ptFN[FNLEN];
  175.     bool    ft;            char ftFN[FNLEN];            tracepos fw;
  176.     bool    rt;            char rtFN[FNLEN];
  177.     bool    ps;            char psFN[FNLEN];
  178.     bool    tfs;        char tfsFN[FNLEN];
  179.     bool    bpd;        char bpdFN[FNLEN];            ushort    neighborhood;
  180.                         char bpdFNPFX[FNLEN];
  181.     bool    prune;        char pruneFN[FNLEN];        tracepos separation;
  182.  
  183.             opts_t(void);
  184.             ~opts_t(void) {}
  185.     void    ClearActions(void);
  186.     int        Digest(int argc, char* argv[]);
  187.     void    MakeTraceFNs(const char* FNprefix);
  188.     void    MakeTracefileFN(const char* FNprefix);
  189.     friend ostream& operator<<(ostream& os, const opts_t& opt);
  190.     };
  191.  
  192. inline
  193. opts_t::opts_t(void):
  194.     pickPeaks(FALSE),
  195.     firstActionSwitch(TRUE),
  196.     format(unknown),
  197.     headers(FALSE),
  198.     verbose(TRUE),
  199.  
  200.     bases((BIT[0]<<A)+(BIT[0]<<C)+(BIT[0]<<G)+(BIT[0]<<T)),
  201.     xform(FALSE),
  202.     scale(FALSE),
  203.     bl(FALSE),
  204.     prune(FALSE),                    separation(SEPARATION_DEFAULT),
  205.     smooth(FALSE),                    smooth_iterations(SMOOTH_ITER_DEFAULT),
  206.     left(0),
  207.     right(0),
  208.     baseline(0),
  209.  
  210.     ind(TRUE),
  211.     d1(TRUE),
  212.     d2(TRUE),
  213.     ts(TRUE),
  214.     pt(TRUE),
  215.     ft(TRUE),                        fw(FLUOR_WIDTH_DEFAULT),
  216.     rt(TRUE),
  217.     tfs(TRUE),
  218.     bp(FALSE),
  219.     fs(FALSE),
  220.     ps(TRUE),
  221.     bpd(FALSE),                        neighborhood(MAX_NBHD),
  222.  
  223.     PMC(PMC_DEFAULT),
  224.     zthresh(Z_THRESH_DEFAULT)
  225.     {
  226.     // ensure that all strings are 0-length
  227.     prefix[0]    =     indFN[0]    =    d1FN[0]        =     d2FN[0]        = \
  228.     tsFN[0]        =    ptFN[0]        =     ftFN[0]        =    rtFN[0]        = \
  229.     tfsFN[0]    =     bpFN[0]        =    fsFN[0]        =    psFN[0]        = \
  230.     xformFN[0]    =    scaleFN[0]    =    bpdFN[0]    =    bpdFNPFX[0]    =
  231.     pruneFN[0]    = NULL;
  232.     }
  233.  
  234. void
  235. opts_t::ClearActions(void)
  236.     // For the first action specified, clear all of the defaults
  237.     {
  238.     if (firstActionSwitch)
  239.         {
  240.         bp = ind = d1 = d2 = fs = ts = pt \
  241.             = ft = rt = ps = tfs = bpd = prune = FALSE;
  242.         firstActionSwitch = FALSE;
  243.         }
  244.     }
  245.  
  246. int
  247. opts_t::Digest(int argc, char** argv)
  248.     {
  249.     uint    argIndex;
  250.     char*    flag;
  251.     bool    hyphenflag = FALSE;            // if we remove a hyphen and get
  252.                                         // something we don't understand, then
  253.                                         // call help
  254.  
  255.     if (argc<2)
  256.         Fatal("bad command line; try autoseq -help");
  257.  
  258.     // scan for version request
  259.     for (argIndex=1;argIndex<=argc-1;argIndex++)
  260.         {
  261.         flag = argv[argIndex];
  262.         if (*flag == HYPHEN)        flag++;        // peel off the hyphen
  263.         if (strcmp(flag,versionsw)==0)
  264.             {
  265.             Banner(cout);
  266.             return DIGEST_NO_PROC;
  267.             }
  268.         }
  269.  
  270.     // scan for help
  271.     for (argIndex=1;argIndex<=argc-1;argIndex++)
  272.         {
  273.         flag = argv[argIndex];
  274.         if (*flag == HYPHEN)        flag++;        // peel off the hyphen
  275.         if (strcmp(flag,helpsw)==0)
  276.             { help(); return DIGEST_NO_PROC; }
  277.         }
  278.  
  279.     // scan for verbose/quiet
  280.     for (argIndex=1;argIndex<=argc-1;argIndex++)
  281.         {
  282.         flag = argv[argIndex];
  283.         if (*flag == HYPHEN)        flag++;        // peel off the hyphen
  284.         if (strcmp(flag,vsw)==0)
  285.             { verbose = TRUE; break; }
  286.         if (strcmp(flag,quietsw)==0)
  287.             { verbose = FALSE; break; }
  288.         }
  289.     
  290.     // if verbose is on, then echo the command line
  291.     //if (verbose)
  292.     //    {
  293.     //    cout << "% ";
  294.     //    for (uint i=0;i<=argc-1;i++)
  295.     //        cout << argv[i] << " ";
  296.     //    cout << "(" << argc-1 << " arguments)" << endl;
  297.     //    }
  298.     
  299.     // Digest flags from left to right until we can't digest any more.  We'll
  300.     // return this value, which hopefully points to the first filename.
  301.     for (argIndex=1; argIndex<=argc-1; argIndex++)
  302.         {
  303.         flag = argv[argIndex];
  304.         hyphenflag = FALSE;
  305.  
  306.         if (*flag == HYPHEN)        { flag++; hyphenflag = TRUE; }
  307.  
  308.         // ignore verbose and verbose
  309.         if ((strcmp(flag,vsw)==0) || (strcmp(flag,quietsw)==0))
  310.             continue;
  311.  
  312.         // ACTIONS
  313.         if (strcmp(flag,indsw)==0)
  314.             { ClearActions(); ind=TRUE; continue; }
  315.         if (strcmp(flag,d1sw)==0)
  316.             { ClearActions(); d1=TRUE; continue; }
  317.         if (strcmp(flag,d2sw)==0)
  318.             { ClearActions(); d2=TRUE; continue; }
  319.         if (strcmp(flag,tssw)==0)
  320.             { ClearActions(); ts=TRUE; continue; }
  321.         if (strcmp(flag,ptsw)==0)
  322.             { ClearActions(); pt=TRUE; continue; }
  323.         if (strcmp(flag,ftsw)==0)
  324.             { ClearActions(); ft=TRUE; continue; }
  325.         if (strcmp(flag,rtsw)==0)
  326.             { ClearActions(); rt=TRUE; continue; }
  327.         if (strcmp(flag,tfssw)==0)
  328.             { ClearActions(); tfs=TRUE; continue; }
  329.         if (strcmp(flag,bpsw)==0)
  330.             { ClearActions(); bp=TRUE; continue; }
  331.         if (strcmp(flag,fssw)==0)
  332.             { ClearActions(); fs=TRUE; continue; }
  333.         if (strcmp(flag,pssw)==0)
  334.             { ClearActions(); ps=TRUE; continue; }
  335.  
  336.         if (strcmp(flag,headsw)==0)
  337.             { headers=TRUE; continue; }
  338.  
  339.         // MODIFIERS
  340.         // These may seg fault if the user omits an argument at the
  341.         // end of the command line
  342.         if (strcmp(flag,pmcsw)==0)
  343.             { PMC = atof(argv[++argIndex]); continue; }
  344.         if (strcmp(flag,basesw)==0)
  345.             {
  346.             bases = 0;
  347.             for (const char *base = argv[++argIndex] ; *base; base++)
  348.                 if (strchr("ACGT",*base)==NULL)
  349.                     Fatal("Invalid base list -- base must be in {A,C,G,T}.");
  350.                 else
  351.                     for (short i=0;i<NUM_BASES;i++)
  352.                         if ((*base == DNA_BASES[i]) && (!(bases>>i & BIT[0])))
  353.                             // base matches && !already selected
  354.                             {bases += (1<<i); break;}
  355.             if (bases == 0)
  356.                 Fatal("Base specification is empty.");
  357.  
  358.             continue;
  359.             }
  360.         if (strcmp(flag,blsw)==0)
  361.             {
  362.             bl = TRUE;
  363.             if (*argv[argIndex+1] != HYPHEN)
  364.                 baseline = atof(argv[++argIndex]);
  365.             continue;
  366.             }
  367.         if (strcmp(flag,bpdsw)==0)
  368.             {
  369.             ClearActions();
  370.             bpd = TRUE;
  371.             if (*argv[argIndex+1] != HYPHEN)    // next argument is a flag
  372.                 {
  373.                 neighborhood = atoi(argv[++argIndex]);
  374.                 if ((neighborhood<1)||(neighborhood>5))
  375.                     Fatal("Neighborhood must be in [1,5].");
  376.                 }
  377.             continue;
  378.             }
  379.         if (strcmp(flag,fmtsw)==0)
  380.             {
  381.             // sequential search through tracefile format abbreviations to find
  382.             // match, or default to unknown format if not found
  383.             ++argIndex;
  384.             short i;
  385.             for (i=0;i<(short)unknown;i++)
  386.                 if (strcmp(TRACEFILE_FORMATS[i].abbr,argv[argIndex])==0)
  387.                     break;
  388.             format = (format_t)i;
  389.             // if the format string is found, then format is set; else = unknown
  390.             continue;
  391.             }
  392.         if (strcmp(flag,fwsw)==0)
  393.             { fw = atoi(argv[++argIndex]); continue; }
  394.         if (strcmp(flag,leftsw)==0)
  395.             {
  396.             if (*argv[argIndex+1]=='p')
  397.                 left = START_AT_PRIMER;
  398.             else
  399.                 left = atoi(argv[argIndex+1]);
  400.             argIndex++;
  401.             continue;
  402.             }
  403.         if (strcmp(flag,outprfxsw)==0)
  404.             { strcpy(prefix,argv[++argIndex]); continue; }
  405.         if (strcmp(flag,prunesw)==0)
  406.             {
  407.             ClearActions();
  408.             prune=TRUE;
  409.             if (*argv[argIndex+1] != HYPHEN)
  410.                 separation = atoi(argv[++argIndex]);
  411.             continue;
  412.             }
  413.         if (strcmp(flag,rightsw)==0)
  414.             { right = atoi(argv[++argIndex]); continue; }
  415.         if (strcmp(flag,scalesw)==0)
  416.             {
  417.             scale=TRUE;
  418.  
  419.             if (*argv[argIndex+1] == HYPHEN)    // next argument is a flag
  420.                 memmove(scales,SCALES,sizeof(SCALES));
  421.             else
  422.                 {
  423.                 strcpy(scaleFN,argv[++argIndex]);
  424.                 ifstream ifs(scaleFN, ios::in);
  425.                 if (ifs.bad())
  426.                     Fatal("couldn't open scale file");
  427.                 ifs >> scales[A] >> scales[C] >> scales[G] >> scales[T];
  428.                 if (ifs.bad())
  429.                     Fatal("error reading scale file");
  430.                 }
  431.             continue;
  432.             }
  433.         if (strcmp(flag,smoothsw)==0)
  434.             {
  435.             smooth = TRUE;
  436.             if ((argIndex<argc-2) && (*argv[argIndex+1] != HYPHEN))
  437.                 smooth_iterations = atoi(argv[++argIndex]);
  438.             continue;
  439.             }
  440.         if (strcmp(flag,xformsw)==0)
  441.             {
  442.             xform=TRUE;
  443.             if (argIndex>argc-2)
  444.                 Fatal("bad command line; try autoseq -help");
  445.             if (*argv[argIndex+1] == HYPHEN)
  446.                 memmove(matrix,ORTHO_MATRIX,sizeof(xform_mtx));
  447.             else
  448.                 {
  449.                 strcpy(xformFN,argv[++argIndex]);
  450.                 cout << "reading " << xformFN << endl;
  451.                 ifstream ifs(xformFN, ios::in);
  452.                 if (ifs.bad())
  453.                     Fatal("couldn't open transformation file");
  454.                 ifs >> matrix[A][A]>>matrix[A][C]>>matrix[A][G]>>matrix[A][T] 
  455.                     >> matrix[C][A]>>matrix[C][C]>>matrix[C][G]>>matrix[C][T] 
  456.                     >> matrix[G][A]>>matrix[G][C]>>matrix[G][G]>>matrix[G][T] 
  457.                     >> matrix[T][A]>>matrix[T][C]>>matrix[T][G]>>matrix[T][T];
  458.                 if (ifs.bad())
  459.                     Fatal("error reading transformation file");
  460.                 }
  461.             continue;
  462.             }
  463.         if (strcmp(flag,zerosw)==0)
  464.             { zthresh = atof(argv[++argIndex]); continue; }
  465.  
  466.         // DEFAULT: unrecognized flags... we're done
  467.         if (hyphenflag)
  468.             // we've stripped off a hyphen, which means the user intended this
  469.             // to be a flag, but we didn't understand it.  Call help and tell
  470.             // our caller to not process.
  471.             {
  472.             char buf[100];
  473.             strcpy(buf,"Unrecognized flag '");
  474.             strcat(buf,flag);
  475.             strcat(buf,"'; try autoseq -help.");
  476.             Fatal(buf);
  477.             }
  478.         break;
  479.         }
  480.  
  481.     // normal exits only
  482.     pickPeaks=(ts || tfs || pt || ft || rt || ps || bpd || prune);
  483.  
  484.     return argIndex;                // filename list starts at argIndex
  485.     }
  486.  
  487. void
  488. opts_t::MakeTraceFNs(const char* FNprefix)
  489.     {
  490.     if (strlen(FNprefix)>FNLEN-SUFXLEN)
  491.       // string too long
  492.       Fatal("Filename is too long to construct destination names.");
  493.  
  494.     char    tempFN[FNLEN];
  495.     strcpy(tempFN,FNprefix);
  496.  
  497.     // Filenames are constructed by replacing the BASE_TOKEN with the base
  498.     // letter.  Thus, every filename which represents a data pertaining to a
  499.     // single base must have a BASE_TOKEN in it.  Check to see if the user
  500.     // provided one; if not, concatenate '.'BASE_TOKEN (ie ".#") to the
  501.     // filename prefix.  The base token will be replaced by the base letter
  502.     // during processing.
  503.     if (NULL == strrchr(tempFN,BASE_TOKEN))
  504.         {
  505.         // BASE_TOKEN not found
  506.         strcat(tempFN,". ");            // leave a space for the BASE_TOKEN
  507.         tempFN[strlen(tempFN)-1]=BASE_TOKEN;
  508.         }
  509.  
  510.     // build filenames for trace-specific data
  511.     if (ind)    { strcpy(indFN,tempFN); }
  512.     if (d1)        { strcpy(d1FN,tempFN); strcat(d1FN,".d1");        }
  513.     if (d2)        { strcpy(d2FN,tempFN); strcat(d2FN,".d2");        }
  514.     if (ts)        { strcpy(tsFN,tempFN); strcat(tsFN,".ts");        }
  515.     if (pt)        { strcpy(ptFN,tempFN); strcat(ptFN,".pt");        }
  516.     if (ft)        { strcpy(ftFN,tempFN); strcat(ftFN,".ft");        }
  517.     if (rt)        { strcpy(rtFN,tempFN); strcat(rtFN,".rt");        }
  518.     }
  519.  
  520. void
  521. opts_t::MakeTracefileFN(const char* FNprefix)
  522.     {
  523.     if (strlen(FNprefix)>FNLEN-SUFXLEN)
  524.       // string too long
  525.       Fatal("Filename is too long to construct destination names.");
  526.  
  527.     char    tempFN[FNLEN];
  528.     strcpy(tempFN,FNprefix);
  529.  
  530.     // remove # if it exists
  531.     char*    start;
  532.     while (start = strrchr(tempFN,BASE_TOKEN))
  533.         // # was found at start... move memory left 1 character
  534.         while (*start = *(start+1))        // shifts memory left until NULL
  535.             start++;
  536.     
  537.     // The filename prefix is now some string, possibly empty, which doesn't
  538.     // contain the BASE_TOKEN.  If the string is not 0-length, then we'll
  539.     // append a period to it to separate the prefix from the suffix to be added
  540.     // subsequently.
  541.     if ((*tempFN != NULL) &&
  542.         (tempFN[strlen(tempFN)-1] != '.') &&
  543.         (tempFN[strlen(tempFN)-1] != '/'))
  544.         strcat(tempFN,".");
  545.  
  546.     // build filenames for files which represent tracefile data
  547.     // (as opposed to trace-specific data)
  548.     if (tfs)    { strcpy(tfsFN,tempFN);        strcat(tfsFN,"tfs");        }
  549.     if (bpd)    {
  550.                 strcpy(bpdFNPFX,tempFN);
  551.                 strcat(bpdFNPFX,"bpd");
  552.                 strcpy(bpdFN,bpdFNPFX);
  553.                 if (neighborhood!=1)
  554.                     sprintf(&bpdFN[strlen(bpdFN)],"[1-%d]",neighborhood);
  555.                 else
  556.                     strcat(bpdFN,"1");
  557.                 }
  558.     if (bp)        { strcpy(bpFN,tempFN);        strcat(bpFN,"bp");            }
  559.     if (fs)        { strcpy(fsFN,tempFN);        strcat(fsFN,"fseq");        }
  560.     if (ps)        { strcpy(psFN,tempFN);        strcat(psFN,"pseq");        }
  561.     if (prune)    { strcpy(pruneFN,tempFN);    strcat(pruneFN,"pruned");    }
  562.     }
  563.  
  564. ostream&
  565. operator<<(ostream& os, const opts_t& opt)
  566.     {
  567.     // general options
  568.         os << "verbose:                  " << YesNo(opt.verbose) << endl;
  569.         os << "headers:                  " << YesNo(opt.headers) << endl;
  570.         os << "file format:              " << FileFormatAbbr(opt.format)
  571.             << " ["    << FileFormatDesc(opt.format) << "]" << endl;
  572.     if (opt.prefix)
  573.         os << "file prefix:              " << opt.prefix << endl;
  574.  
  575.     // baseline
  576.     if (opt.bl)
  577.         {
  578.         os << "translate traces by:      ";
  579.         if (opt.baseline==0)
  580.                                         os << "min value";
  581.         else
  582.                                         os << -opt.baseline;
  583.         os << endl;
  584.         }
  585.  
  586.     // scaling
  587.     if (opt.scale)
  588.         {
  589.         os << "scale A,C,G,T traces by";
  590.         if (NULL != opt.scaleFN)
  591.             os << " '" << opt.scaleFN << "'";
  592.         os << ":" << endl;
  593.         os     << "  A: " << opt.scales[A] << endl
  594.             << "  C: " << opt.scales[C] << endl
  595.             << "  G: " << opt.scales[G] << endl
  596.             << "  T: " << opt.scales[T] << endl;
  597.         }
  598.  
  599.     // smoothing
  600.     if (opt.smooth)
  601.         os << "smooth traces:            "
  602.             << opt.smooth_iterations
  603.             << " time" << Plural(opt.smooth_iterations) << endl;
  604.  
  605.     // transform matrix
  606.     if (opt.xform)
  607.         {
  608.         os << "transform trace vector with ";
  609.         if (NULL == *opt.xformFN)
  610.             os << "standard matrix." << endl;
  611.         else
  612.             os << "'" << opt.xformFN << "'." << endl;
  613.         os    << "          A        C        G        T" << endl;
  614.         for (uint i=A;i<NUM_BASES;i++)
  615.             {
  616.             cout << DNA_BASES[i] << "   ";
  617.             for (uint j=A;j<NUM_BASES;j++)
  618.                 {
  619.                 os.setf(ios::right,ios::adjustfield);
  620.                 os.width(7);
  621.                 os.precision(3);
  622.                 os << opt.matrix[i][j] << "  ";
  623.                 }
  624.             cout << endl;
  625.             }
  626.         }
  627.  
  628.     // actions to take, and specification of bases and trace indices
  629.     if (opt.ind || opt.d1 || opt.d2 || opt.ts || opt.pt || opt.ft || opt.rt)
  630.                     // only write base specification if it'll be used
  631.                     os << "for bases "
  632.                         << (opt.bases>>A & BIT[0] ? "A" : "")
  633.                         << (opt.bases>>C & BIT[0] ? "C" : "")
  634.                         << (opt.bases>>G & BIT[0] ? "G" : "")
  635.                         << (opt.bases>>T & BIT[0] ? "T" : "")
  636.                         << " ";
  637.     if ((0 == opt.left) && (0 == opt.right))
  638.                     os << "over entire trace range";
  639.     else
  640.         {
  641.         // print range, with handling for special cases
  642.         os    << "in trace range [";
  643.         if (opt.left==START_AT_PRIMER)
  644.             os << "primer";
  645.         else
  646.             os << opt.left;
  647.         os << ",";
  648.         if (opt.right<1)
  649.             os << -opt.right << " from end";
  650.         else
  651.             os << opt.right;
  652.         os << "]";
  653.         }
  654.     os << ", do the following:" << endl;
  655.  
  656.     if (opt.ind) os<<"    individual traces to:     '" << opt.indFN <<"'"<<endl;
  657.     if (opt.d1)  os<<"    1st derivative traces to: '" << opt.d1FN  <<"'"<<endl;
  658.     if (opt.d2)  os<<"    2nd derivative traces to: '" << opt.d2FN  <<"'"<<endl;
  659.     if (opt.ts)  os<<"    trace summaries to:       '" << opt.tsFN  <<"'"<<endl;
  660.     if (opt.pt)  os<<"    peak traces to:           '" << opt.ptFN  <<"'"<<endl;
  661.     if (opt.ft)  os<<"    fluorescence traces to:   '" << opt.ftFN  <<"'"<<endl;
  662.     if (opt.rt)  os<<"    residual traces to:       '" << opt.rtFN  <<"'"<<endl;
  663.     if (opt.bp)  os<<"    base positions to:        '" << opt.bpFN  <<"'"<<endl;
  664.     if (opt.fs)  os<<"    file sequence to:         '" << opt.fsFN  <<"'"<<endl;
  665.     if (opt.ps)  os<<"    predicted sequence to:    '" << opt.psFN  <<"'"<<endl;
  666.     if (opt.tfs) os<<"    tracefile summary to:     '" << opt.tfsFN <<"'"<<endl;
  667.     if (opt.bpd) os<<"    base pos. deltas ("
  668.                    << opt.neighborhood << ") to:  '"   << opt.bpdFN <<"'"<<endl;
  669.  
  670.     if (opt.pickPeaks)
  671.         // only write these if they'll be used
  672.         {
  673.         if (opt.prune) os<<"    pruned peak list ("
  674.                    << opt.separation   << ") to:  '" << opt.pruneFN <<"'"<<endl;
  675.         os << "Peak height threshold =   " << opt.PMC
  676.             << " * Mean Trace Value" << endl;
  677.         os << "1d zero threshold =       " << opt.zthresh << endl;
  678.         os << "Gaussian model width =    " << 2 * opt.fw + 1
  679.             << " (= center pt + 2 * " << opt.fw
  680.             << " pts each side of center)" << endl;
  681.                     }
  682.     return os;
  683.     }
  684.  
  685. inline
  686. char*
  687. Substitute(const char* src, char* dest, uint base)
  688.     // substitute base letter for place holder in filename
  689.     // I assume that BASE_TOKEN is in the string, which is guaranteed for the
  690.     // use in this source... and guaranteed to seg fault if this assumption is
  691.     // false.
  692.     {
  693.     strcpy(dest,src);
  694.     *strrchr(dest,BASE_TOKEN)=DNA_BASES[base];
  695.     return dest;
  696.     }
  697.  
  698.  
  699. #ifdef WIN_MAC
  700. //#include <console.h>
  701. extern "C" int ccommand(char **argv[]);
  702.  
  703. // use with ChildApp.c
  704. extern "C" int RealMain( int argc, char *argv[])
  705.  
  706. #else
  707.  
  708. int
  709. main(int argc, char* argv[])
  710. #endif
  711.     {
  712.     opts_t        options;
  713.     int            argIndex;                    // the next arg to read
  714.     char        outpath[FNLEN];
  715.  
  716. #ifdef WIN_MAC
  717.     //argc= ccommand( &argv);
  718. #endif
  719.  
  720.     argIndex = options.Digest(argc, argv);
  721.  
  722.     if (argIndex == DIGEST_NO_PROC)            // we've been told to not process
  723.         return 0;
  724.  
  725.     if (options.verbose)
  726.         Banner(cout);
  727.     
  728.     while (argIndex<argc)
  729.         {
  730.         if (options.prefix[0] == NULL)
  731.             {
  732.             // no output specified... prefix is just input filename
  733.             // strcpy(outpath,argv[argIndex]);
  734.  
  735.             // no output specified... prefix is empty
  736.             // (actually, using prefix = # means base, and will be trimmed to
  737.             // empty for output which isn't base-specific)
  738.             outpath[0]=BASE_TOKEN;
  739.             outpath[1]=NULL;
  740.             }
  741.         else
  742.             {
  743.             // output specified... if it ends with '/', then it's a directory
  744.             // prefix and we need to append the filename to it; otherwise, it's
  745.             // a full file prefix specifier
  746.             strcpy(outpath,options.prefix);
  747.             if (outpath[strlen(outpath)-1]=='/')
  748.                 {
  749.                 //strcat(outpath,argv[argIndex]);
  750.                 strcat(outpath,".");
  751.                 outpath[strlen(outpath)-1]=BASE_TOKEN;
  752.                 }
  753.             }
  754.         // by now, outpath is a full file pathname
  755.         // chop off the .s1 suffix if it exists
  756. // dgg mod >> need to generalize this beyond ".s1" or leave s1 on...
  757. #if 1
  758.  
  759. #else
  760.         if ((strlen(outpath)>=3) &&
  761.             (strcmp(&outpath[strlen(outpath)-3],".s1")==0))
  762.             outpath[strlen(outpath)-3]=NULL;
  763. #endif
  764.  
  765.         options.MakeTraceFNs(outpath);
  766.         options.MakeTracefileFN(outpath);
  767.     
  768.         if (options.verbose)
  769.             cout << setfill('=') << setw(80) << "" << setfill(' ') << endl
  770.                 << options;
  771.     
  772.         if (process(argv[argIndex],options))
  773.             Fatal("Error occurred during processing");
  774.  
  775.         argIndex++;
  776.         }
  777.  
  778.     if (options.verbose)
  779.         cout << setfill('=') << setw(80) << "" << setfill(' ') << endl;
  780.  
  781.     return 0;
  782.     }
  783.  
  784. int
  785. process(const char *FN,    opts_t& options)
  786.     {
  787.     const bool&    verbose = options.verbose;
  788.     const bool&    headers = options.headers;
  789.     const short& bases = options.bases;
  790.     format_t    fmt = options.format;
  791.     uint        traceIndex;
  792.     char        tempFN[FNLEN];
  793.     //ofstream*    os;
  794.     CTraceFile::error_t CTFError = CTraceFile::noError;
  795.     CTError        CTError = noError;
  796.     //CTraceFile     tracefile(FN,fmt);        // declared when file is read
  797.  
  798.     //
  799.     // read file
  800.     //
  801.     if (fmt == unknown)
  802.         fmt = FileFormat(FN);
  803.  
  804.     if (verbose)
  805.         cout << setfill('-') << setw(80) << "" << setfill(' ') << endl
  806.             <<"Reading '"<<FN<<"' ("<<FileFormatAbbr(fmt)<< ")..." << flush;
  807.  
  808.     CTraceFile     tracefile(FN,fmt);
  809.     CTFError = tracefile.Error();
  810.     if (verbose)
  811.         cout << "(" << (int)CTFError << ")." << endl;
  812.  
  813.     // assume file doesn't exist if fmt is FFerr
  814.     if (CTraceFile::noError != CTFError)
  815.         {
  816.         char buf[200];
  817.         strcpy(buf,"error while reading file '");
  818.         strcat(buf,FN);
  819.         strcat(buf,"'; trying to continue.");
  820.         Warning(buf);
  821.         return 0; // continue
  822.         }
  823.  
  824.     //
  825.     // set cutoffs
  826.     //
  827.     if (options.left !=0)
  828.         {
  829.         tracepos lc = (    options.left==START_AT_PRIMER ?
  830.                         tracefile.PrimerPosition() :
  831.                         options.left);
  832.         if (lc<0)
  833.             Fatal("Bad left cutoff specification; check trace limits.");
  834.         if (lc>tracefile.NumPoints()-2)
  835.             Fatal("Left cutoff would leave empty trace window.");
  836.         tracefile.LeftCutoff(lc);
  837.         }
  838.     if (options.right!=0)
  839.         {
  840.         // remember... options.right<0, therefore we're shifting it /left/
  841.         tracepos rc = (    options.right>0 ?
  842.                         options.right :
  843.                         tracefile.NumPoints()+options.right);
  844.         if (rc<tracefile.LeftCutoff())
  845.                 Fatal("Attempt to right cutoff < left cutoff.");
  846.         if (rc>tracefile.RightCutoff())
  847.                 Fatal("Attempt to set right cutoff > size of trace.");
  848.     
  849.         tracefile.RightCutoff(rc);
  850.         }
  851.  
  852.  
  853.     if (verbose)
  854.         cout << tracefile
  855.             << setfill('-') << setw(80) << "" << setfill(' ') << endl;
  856.  
  857.     //
  858.     // smooth
  859.     //
  860.     if (options.smooth)
  861.         {
  862.         if (verbose)
  863.             cout << "smoothing ("
  864.                 << options.smooth_iterations
  865.                 << " passes):\t";
  866.  
  867.         ushort iter = options.smooth_iterations;
  868.         while(iter--)
  869.             for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  870.                 tracefile[traceIndex]->Smooth();
  871.  
  872.         tracefile.CalculateStats();
  873.  
  874.         if (verbose)
  875.             cout << "done." << endl;
  876.         }
  877.  
  878.     //
  879.     // translate
  880.     //
  881.     if (options.bl)
  882.         {
  883.         if (verbose)
  884.             cout << "translating: (";
  885.  
  886.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  887.             {
  888.             double bl = options.baseline;
  889.             if (0 == bl)
  890.                 bl = tracefile[traceIndex]->Min();
  891.  
  892.             if (verbose)
  893.                 cout <<DNA_BASES[traceIndex]<<":"<<-bl;
  894.             if (verbose && (traceIndex!=T)) cout <<"; "<<flush;
  895.  
  896.             tracefile[traceIndex]->Translate(-bl);
  897.             }
  898.  
  899.         if (verbose)
  900.             cout << ")" << endl;
  901.         }
  902.  
  903.     //
  904.     // scaling
  905.     //
  906.     if (options.scale)
  907.         {
  908.         if (verbose)
  909.             cout << "scaling traces:\t";
  910.  
  911.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  912.             {
  913.             if (verbose)
  914.                 cout << 1/options.scales[traceIndex] << " " <<flush;
  915.  
  916.             tracefile[traceIndex]->Scale(1/options.scales[traceIndex]);
  917.             }
  918.  
  919.         if (verbose)
  920.             cout << endl;
  921.         }
  922.  
  923.     //
  924.     // transform
  925.     //
  926.     if (options.xform)
  927.         {
  928.         if (verbose)
  929.             cout << "transforming traces..." << flush;
  930.  
  931.         CTFError = tracefile.Transform(options.matrix);
  932.  
  933.         if (verbose)
  934.             cout << "(" << (int)CTFError << "); " << flush;
  935.  
  936.         if (CTraceFile::noError != CTFError)
  937.             Fatal("Transformation failed...exiting");
  938.  
  939.         if (verbose)
  940.             cout << "translating: (";
  941.  
  942.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  943.             {
  944.             double bl = tracefile[traceIndex]->Min();
  945.             if (verbose)    cout <<DNA_BASES[traceIndex]<<":"<<-bl;
  946.             if (verbose && (traceIndex!=T)) cout <<"; "<<flush;
  947.             tracefile[traceIndex]->Translate(-bl);
  948.             }
  949.  
  950.         if (verbose)
  951.             cout << ")" << endl;
  952.         }
  953.  
  954.     //
  955.     // write individual traces
  956.     //
  957.     if (options.ind)
  958.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  959.             {
  960.             if (!(bases>>traceIndex & BIT[0]))    continue;
  961.  
  962.             Substitute(options.indFN,tempFN,traceIndex);
  963.             //os = new ofstream(tempFN,ios::out);
  964.             ofstream os(tempFN);  
  965.  
  966.             if (verbose)
  967.                 cout << "Writing individual trace " << DNA_BASES[traceIndex]
  968.                     << " to '" <<tempFN<<"'..."<<flush;
  969.  
  970.             if (headers)
  971.                 //*os << Banner
  972.                 os << Banner
  973.                     << COMMENT_PREFIX << DNA_BASES[traceIndex] << " Trace of "
  974.                     << "'"<<FN<<"'" << endl
  975.                     << *tracefile[traceIndex];
  976.  
  977.             CTError = tracefile[traceIndex]->WriteAsText(os); //*os 
  978.  
  979.             if (verbose)
  980.                 cout << "(" << (int)CTError << ")." << endl;
  981.  
  982.             //delete os;
  983.  
  984.             if (noError != CTError)
  985.                 Fatal("Error writing individual trace.");
  986.             }
  987.  
  988.     //
  989.     // write derivative traces
  990.     //
  991.     if (options.d1 || options.d2)
  992.         {
  993.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  994.             {
  995.             if (!(bases>>traceIndex & BIT[0]))    continue; // skip this base
  996.  
  997.             CTrace<trace_t>&    trace = *tracefile[traceIndex];
  998.             CTrace<double>*    d1trace=NULL;
  999.             CTrace<double>*    d2trace=NULL;
  1000.  
  1001.             // make the 1d trace (need for the d2 even if we don't write d1)
  1002.             if (verbose)
  1003.                 cout << "Trace " << DNA_BASES[traceIndex]
  1004.                     << ": making 1st derivative..." << flush;
  1005.  
  1006.             CTError = trace.Derivative(&d1trace);
  1007.  
  1008.             if (verbose)
  1009.                 cout << "(" << (int)CTError << "); " << flush;
  1010.  
  1011.             if (noError != CTError)
  1012.                 {
  1013.                 if (verbose) cout << endl;
  1014.                 Fatal("Error making derivative.");
  1015.                 }
  1016.  
  1017.             // write the d1 trace if we're supposed to (we may have created
  1018.             // d1 just to make d2)
  1019.             if (options.d1)
  1020.                 {
  1021.                 Substitute(options.d1FN,tempFN,traceIndex);
  1022.  
  1023.                 if (verbose)
  1024.                      cout << "writing '" << tempFN << "'..." << flush;
  1025.  
  1026.                 //os = new ofstream(tempFN,ios::out);
  1027.                 ofstream os(tempFN);
  1028.                 if (headers)
  1029.                     os << Banner
  1030.                         << COMMENT_PREFIX << "1st derivative of "
  1031.                         << "'"<<FN<<"'" << endl
  1032.                         << *d1trace;
  1033.                 CTError = d1trace->WriteAsText(os);
  1034.  
  1035.                 if (verbose)
  1036.                     cout << "(" << (int)CTError << ")." << endl;
  1037.  
  1038.                 if (noError != CTError)
  1039.                     Fatal("Error writing derivative.");
  1040.                 } // options.d1
  1041.             else
  1042.                 if (verbose)
  1043.                     cout << "not written." << endl;
  1044.  
  1045.             // make and write the d2 trace
  1046.             if (options.d2)
  1047.                 {
  1048.                 if (verbose)
  1049.                     cout << "Trace " << DNA_BASES[traceIndex]
  1050.                         << ": making 2nd derivative..." << flush;
  1051.  
  1052.                 CTError = d1trace->Derivative(&d2trace);
  1053.  
  1054.                 if (verbose)
  1055.                     cout << "(" << (int)CTError << "); " << flush;
  1056.  
  1057.                 if (noError != CTError)
  1058.                     {
  1059.                     if (verbose) cout << endl;
  1060.                     Fatal("Error making derivative.");
  1061.                     }
  1062.                 Substitute(options.d2FN,tempFN,traceIndex);
  1063.  
  1064.                 if (verbose)
  1065.                     cout << "writing '" << tempFN << "'..." << flush;
  1066.  
  1067.                 ofstream os(tempFN);
  1068.                 if (headers)
  1069.                     os << Banner
  1070.                         << COMMENT_PREFIX << "2nd derivative of "
  1071.                         << "'"<<FN<<"'"
  1072.                         << *d2trace << endl;
  1073.                 CTError = d2trace->WriteAsText(os);
  1074.  
  1075.                 if (verbose)
  1076.                     cout << "(" << (int)CTError << ")." << endl;
  1077.  
  1078.                 if (noError != CTError)
  1079.                     Fatal("Error writing derivative.");
  1080.                 delete d2trace;
  1081.                 } // options.d2
  1082.             delete d1trace;
  1083.             } // for loop through traces
  1084.         }
  1085.         // rof (trace loop)
  1086.     // fi (d1 || d2)
  1087.  
  1088.     //
  1089.     // pick peaks
  1090.     //
  1091.     if (options.pickPeaks)
  1092.         {
  1093.         // set baseline values
  1094.         if (verbose)
  1095.             cout << "Setting baseline to " << options.baseline
  1096.                 << " for all traces..." << flush;
  1097.  
  1098.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  1099.             tracefile[traceIndex]->Baseline(options.baseline);
  1100.  
  1101.         if (verbose)
  1102.             cout << "done." << endl;
  1103.  
  1104.         if (verbose)
  1105.             cout << "Picking peaks for all traces..." << flush;
  1106.  
  1107.         CTFError = tracefile.PickPeaks(options.PMC,options.zthresh);
  1108.  
  1109.         if (verbose)
  1110.             cout << "(" << (int)CTFError << ").  ["
  1111.                 << tracefile.Peaks().Size() << " peaks picked]"
  1112.                 << endl;
  1113.  
  1114.         if (CTFError != CTraceFile::noError)
  1115.             Fatal("Error during peakpicking.");
  1116.  
  1117.  
  1118.         //
  1119.         // adjust all peak records to record positions relative to leftCutoff
  1120.         //
  1121.         if (tracefile.LeftCutoff() != 0)
  1122.             {
  1123.             // translate the tracefile's peak records
  1124.             tracefile.Peaks().Offset(-tracefile.LeftCutoff());
  1125.  
  1126.             // and the peak records for the individual traces
  1127.             for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  1128.                 tracefile[traceIndex]->Peaks().Offset(-tracefile.LeftCutoff());
  1129.             }
  1130.  
  1131.  
  1132.         //
  1133.         // prune peaks
  1134.         //
  1135.         if (options.prune)
  1136.             {
  1137.             ofstream os(options.pruneFN);
  1138.  
  1139.             if (verbose)
  1140.                 cout << "Pruning peak list (sep=" << options.separation<< ")\n"
  1141.                      << "   writing to '" << options.pruneFN << "'..." <<flush;
  1142.  
  1143.             ulong    szBefore = tracefile.Peaks().Size();
  1144.  
  1145.             if (headers)
  1146.                 os    << "minimum pairwise separation = " << options.separation
  1147.                     << "; deleted peaks marked with *; output is 100 cols wide"
  1148.                     << endl
  1149.                     << PeakRecHeader;
  1150.  
  1151.             tracefile.PrunePeaks(options.separation, &os);
  1152.  
  1153.             if (verbose)
  1154.                 cout << szBefore << "-"<< tracefile.Peaks().Size()
  1155.                      << "=" << szBefore-tracefile.Peaks().Size()
  1156.                      << " peaks removed." << endl;
  1157.             }
  1158.  
  1159.         //
  1160.         // write individual peak traces
  1161.         //   (do this before shifting peak records in the next step)
  1162.         if (options.pt)
  1163.             for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  1164.                 {
  1165.                 if (!(bases>>traceIndex & BIT[0]))    continue;
  1166.                 CTrace<trace_t>&    trace = *tracefile[traceIndex];
  1167.  
  1168.                 Substitute(options.ptFN,tempFN,traceIndex);
  1169.  
  1170.                 if (verbose)
  1171.                     cout << "Trace " << DNA_BASES[traceIndex]
  1172.                         << " peak trace: computing..." << flush;
  1173.  
  1174. //                trace.Peaks().ComputePTrace(trace.Size());
  1175.                 trace.Peaks().ComputePTrace(trace);
  1176.  
  1177.                 if (verbose)
  1178.                     cout << "writing to '" << tempFN << "'..." << flush;
  1179.  
  1180.                 ofstream os(tempFN);
  1181.                 if (headers)
  1182.                     os << Banner
  1183.                         << COMMENT_PREFIX << "Peak trace "
  1184.                         << DNA_BASES[traceIndex] << " of file "
  1185.                         << "'" << FN << "'" << endl;
  1186.  
  1187.                 CTError = trace.Peaks().PTrace->WriteAsText(os);
  1188.  
  1189.                 if (verbose)
  1190.                     cout << "(" << (int)CTError << ")." << endl;
  1191.                 }
  1192.         }
  1193.  
  1194.     //
  1195.     // write individual trace summaries
  1196.     //
  1197.     if (options.ts)
  1198.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  1199.             {
  1200.             if (!(bases>>traceIndex & BIT[0]))    continue;
  1201.             CTrace<trace_t>&    trace = *tracefile[traceIndex];
  1202.     
  1203.             Substitute(options.tsFN,tempFN,traceIndex);
  1204.  
  1205.             if (verbose)
  1206.                 cout << "Writing trace " << DNA_BASES[traceIndex]
  1207.                     << " summary (" << trace.Peaks().Size()
  1208.                     << " peaks) to '" << tempFN << "'..." << flush;
  1209.  
  1210.             ofstream os(tempFN);
  1211.             if (headers)
  1212.                 os << Banner
  1213.                     << COMMENT_PREFIX << "Summary of trace "
  1214.                     << DNA_BASES[traceIndex] << " of file " << "'"<<FN<<"'"
  1215.                     << endl;
  1216.  
  1217.             os  
  1218. #ifndef WIN_MAC
  1219.                 << trace  
  1220. #endif
  1221.                 << trace.Peaks();
  1222.                 
  1223.             if (verbose)
  1224.                 cout << "done." << endl;
  1225.             }
  1226.  
  1227.     //
  1228.     // write individual fluorescence traces
  1229.     // (already computed from PickPeaks)
  1230.     //
  1231.     if (options.ft)
  1232.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  1233.             {
  1234.             if (!(bases>>traceIndex & BIT[0]))    continue;
  1235.             CTrace<trace_t>&    trace = *tracefile[traceIndex];
  1236.  
  1237.             Substitute(options.ftFN,tempFN,traceIndex);
  1238.  
  1239.             if (verbose)
  1240.                 cout << "Trace " << DNA_BASES[traceIndex]
  1241.                     << " expected fluorescence: writing to '"
  1242.                     << tempFN << "'..." << flush;
  1243.  
  1244.             ofstream os(tempFN);
  1245.             if (headers)
  1246.                 os << Banner
  1247.                     << COMMENT_PREFIX << "Fluorescence trace "
  1248.                     << DNA_BASES[traceIndex] << " of file " << "'"<<FN<<"'"
  1249.                     << endl;
  1250.             CTError = trace.Peaks().FTrace->WriteAsText(os);
  1251.  
  1252.             if (verbose)
  1253.                 cout << "(" << (int)CTError << ")." << endl;
  1254.             }
  1255.  
  1256.     //
  1257.     // write individual residual traces
  1258.     //
  1259.     if (options.rt)
  1260.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  1261.             {
  1262.             if (!(bases>>traceIndex & BIT[0]))    continue;
  1263.             CTrace<trace_t>&    trace = *tracefile[traceIndex];
  1264.  
  1265.             Substitute(options.rtFN,tempFN,traceIndex);
  1266.  
  1267.             if (verbose)
  1268.                 cout << "Trace " << DNA_BASES[traceIndex]
  1269.                     << " residual: writing to '"
  1270.                     << tempFN << "'..." << flush;
  1271.  
  1272.             ofstream os(tempFN);
  1273.             if (headers)
  1274.                 os << Banner
  1275.                     << COMMENT_PREFIX << "Residual trace "
  1276.                     << DNA_BASES[traceIndex] << " of file '" << FN << "'"<<endl;
  1277.             CTError = trace.Peaks().RTrace->WriteAsText(os);
  1278.  
  1279.             if (verbose)
  1280.                 cout << "(" << (int)CTError << ")." << endl;
  1281.             }
  1282.  
  1283.     //
  1284.     // write tracefile summary
  1285.     //
  1286.     if (options.tfs)
  1287.         {
  1288.         if (verbose)
  1289.             cout << "Writing tracefile summary (" << tracefile.NumPeaks()
  1290.                 << " peaks) to '" << options.tfsFN << "'..." << flush;
  1291.  
  1292.         ofstream os(options.tfsFN);
  1293.         if (headers)
  1294.             os << Banner
  1295.                 << COMMENT_PREFIX << "Summary of file '" << FN <<"'"<< endl;
  1296.         os << tracefile;
  1297.  
  1298.         if (verbose)
  1299.             cout << "done." << endl;
  1300.         }
  1301.  
  1302.     //
  1303.     // write base position deltas
  1304.     //
  1305.     if (options.bpd)
  1306.         for (ushort nbhd=1; nbhd<=options.neighborhood; nbhd++)
  1307.             {
  1308.             // bpd filename generated by appending the nbhd to bpd
  1309.             strcpy(tempFN,options.bpdFNPFX);
  1310.             sprintf(&tempFN[strlen(tempFN)],"%d",nbhd);
  1311.  
  1312.             if (verbose)
  1313.                 cout << "Writing base position deltas for nbhd=" << nbhd
  1314.                      << " to '" << tempFN << "'..." << flush;
  1315.  
  1316.             ofstream os(tempFN);
  1317.             if (headers)
  1318.                 os << Banner
  1319.                     << COMMENT_PREFIX << "Base position deltas of file '"
  1320.                     << FN << "' for nbhd=" << nbhd << endl;
  1321.  
  1322.             tracefile.Peaks().WriteDeltas(os,nbhd,headers);
  1323.  
  1324.             if (verbose)
  1325.                 cout << "done." << endl;
  1326.             }
  1327.  
  1328.  
  1329.     //
  1330.     // write file sequence
  1331.     //
  1332.     if (options.fs)
  1333.         {
  1334.         if (verbose)
  1335.             cout << "Writing file sequence (" << strlen(tracefile.Sequence())
  1336.                 << " bases) to '" << options.fsFN << "'..." << flush;
  1337.  
  1338.         ofstream os(options.fsFN);
  1339.         if (headers)
  1340.             os << Banner
  1341.                 << COMMENT_PREFIX << "Sequence of file '" << FN << "'" << endl;
  1342.         os << tracefile.Sequence();
  1343.  
  1344.         if (verbose)
  1345.             cout << "done." << endl;
  1346.         }
  1347.  
  1348.     //
  1349.     // write predicted sequence
  1350.     //
  1351.     if (options.ps)
  1352.         {
  1353.         CSequence<PeakRec>&    peaks = tracefile.Peaks();
  1354.  
  1355.         if (verbose)
  1356.             cout << "Writing predicted sequence (" << tracefile.Peaks().Size()
  1357.                 << " bases) to '" << options.psFN << "'..." << flush;
  1358.  
  1359.         ofstream os(options.psFN);
  1360.         if (headers)
  1361.             os << Banner
  1362.                 << COMMENT_PREFIX << "Predicted sequence of tracefile '" << FN
  1363.                     << "'" << endl;
  1364.         char* seq=tracefile.Peaks().Sequence();
  1365.         os << seq;
  1366.         delete seq;
  1367.  
  1368.         if (verbose)
  1369.             cout << "done." << endl;
  1370.         }
  1371.  
  1372.     //
  1373.     // write base positions
  1374.     //
  1375.     if (options.bp)
  1376.         {
  1377.         if (verbose)
  1378.             cout << "Writing base positions (" << strlen(tracefile.Sequence())
  1379.                 << " bases) to '" << options.bpFN << "..." << flush;
  1380.  
  1381.         ofstream os(options.bpFN);
  1382.         if (headers)
  1383.             os << Banner
  1384.                 <<COMMENT_PREFIX<<"DNA base positions of file '"<<FN<<"'"<<endl;
  1385.         for (uint i=0;i<strlen(tracefile.Sequence());i++)
  1386.             os << tracefile.BasePositions()[i] << "\t"
  1387.                 << tracefile.Sequence()[i] << endl;
  1388.  
  1389.         if (verbose)
  1390.             cout << "done." << endl;
  1391.         }
  1392.  
  1393.     return 0;
  1394.     }
  1395.  
  1396. void
  1397. help(void)
  1398.     {
  1399. const int    FLAG_FLD_WIDTH    =4;
  1400. #define        INDENT            "    "        // FLG_FLD_WIDTH
  1401.  
  1402. // customizable delimiters to facilitate exporting to man page, FrameMaker table
  1403. // in order to make these concatenatable, they're defines rather than const
  1404. // the output format is:
  1405. //        [FP <flag> FS AP [<arglist> AS] DP <first line> [LS <lines>] ER]
  1406. // [] = repeated/optional element
  1407.  
  1408.  
  1409. #if FALSE
  1410.     // tab-delimited fields, return delimited records for flags
  1411.     #define FP    << ""
  1412.     #define FS    ""
  1413.     #define    AP    "\t"
  1414.     #define AS    ""
  1415.     #define DP    "\t"
  1416.     #define LS    " "
  1417.     #define ER    "\n"
  1418. #elif FALSE
  1419.     // man page output
  1420.     #define FP    << setw(FLAG_FLD_WIDTH)    << ".B \-"
  1421.     #define FS    " "
  1422.     #define    AP    ""
  1423.     #define AS    "\n" INDENT    FS
  1424.     #define DP    ""
  1425.     #define LS    AS
  1426.     #define ER    "\n"
  1427. #else
  1428.     // standard output
  1429.     #define FP    << setw(FLAG_FLD_WIDTH)    // Flag Prefix; ensures std flag width
  1430.     #define FS    " "                        // Flag Suffix; at least one space after
  1431.     #define    AP    ""                        // Arg Prefix 
  1432.     #define AS    "\n" INDENT    FS            // Arg Suffix
  1433.     #define DP    ""                        // Description Prefix
  1434.     #define LS    AS                        // Line Separator
  1435.     #define ER    "\n"                    // End Record
  1436. #endif                                    // end delimiter definitions
  1437.  
  1438.  
  1439. cout.setf(ios::left);
  1440.  
  1441. cout
  1442. << "autoseq (version " << VERSION << ") help\n"
  1443. << "-------------------------------\n"
  1444. << "  autoseq extracts a variety of information from SCF and ABI files and\n"
  1445. << "  writes these components to files.  It acts as an interface to CTrace,\n"
  1446. << "  CTraceFile, & CPeakList to perform a operations on the data, and the\n"
  1447. << "  results of these operations are also directed to files.\n\n";
  1448.  
  1449. cout
  1450. << "command line structure:\n"
  1451. << "% autoseq {flag_list} {filename}\n"
  1452. << "flag_list ::= {'-'} {flag} {flag_list}\n"
  1453. << "flag ::= 1 of the following:\n";
  1454.  
  1455. cout
  1456. FP << basesw << FS AP "[A,C,G,T]0+" AS
  1457. DP "Specifies to which bases individual actions should be performed." 
  1458. LS "More than one base may be specified. The default is -bACGT (all bases)." 
  1459. << ER
  1460.  
  1461. FP << blsw << FS AP "[short bl=<trace minimum>]" AS
  1462. DP "Set the baseline of the selected traces to bl. If bl is omitted then each" 
  1463. LS "trace is translated downward by the the minimum value for that trace."
  1464. << ER
  1465.  
  1466. FP << bpdsw << FS AP "[+integer maxnbhd=" << MAX_NBHD << "]" AS
  1467. DP "For every subsequence of length nbhd (1<=nbhd<=maxnbhd) in the predicted"
  1468. LS "sequence, write the index of the 3' most peak center for that subsequence,"
  1469. LS "the delta for the subsequence, and the subsequence itself.  The delta"
  1470. LS "between bases m and n is defined to be the distance between the centers of"
  1471. LS "the peaks corresponding to bases m and n."
  1472. << ER
  1473.  
  1474. FP << bpsw << FS AP
  1475. DP "Output the base-to-trace-position mapping stored in file [prefix.bp]" 
  1476. << ER
  1477.  
  1478. FP << d1sw << FS AP
  1479. DP "Compute 1st derivative trace [prefix.#.d1]" 
  1480. << ER
  1481.  
  1482. FP << d2sw << FS AP
  1483. DP "Compute 2nd derivative trace [prefix.#.d2]" 
  1484. << ER
  1485.  
  1486. FP << fmtsw << FS AP "[{ABI0,ABI1,ABI,SCF} fmt=<best guess>]" AS
  1487. DP "Read in specified format; guess if not specified." 
  1488. LS INDENT << setw(5) << FileFormatAbbr(ABI0) << " - " << FileFormatDesc(ABI0) <<
  1489. LS INDENT << setw(5) << FileFormatAbbr(ABI1) << " - " << FileFormatDesc(ABI1) <<
  1490. LS INDENT << setw(5) << FileFormatAbbr(ABI)  << " - " << FileFormatDesc(ABI)  <<
  1491. LS INDENT << setw(5) << FileFormatAbbr(SCF)  << " - " << FileFormatDesc(SCF)
  1492. << ER
  1493.  
  1494. FP << fssw << FS AP
  1495. DP "Write the sequence predicted stored in the file. [prefix.fseq]"
  1496. << ER
  1497.  
  1498. FP << ftsw << FS AP
  1499. DP "Writes the trace of expected fluorescence from the list of selected peaks."
  1500. LS "Expected fluorescence is calculated by fitting a Gaussian at the peak"
  1501. LS "center using the peak height & width; see "<<fwsw<<". [prefix.#.ft]" 
  1502. << ER
  1503.  
  1504. FP << fwsw << FS AP "[+short fw=" << FLUOR_WIDTH_DEFAULT << "]" AS
  1505. DP "Specifies the breadth of the Gaussian used to model the fluorescence of" 
  1506. LS "each peak (ie. Gaussian limits are [peakcenter - fw , peakcenter + fw])." 
  1507. << ER
  1508.  
  1509. FP << headsw << FS AP
  1510. DP "Include descriptive headers in output files" 
  1511. << ER
  1512.  
  1513. FP << helpsw << FS AP
  1514. DP "This help display" 
  1515. << ER
  1516.  
  1517. FP << indsw << FS AP
  1518. DP "Write the individual traces after translation ("<<blsw<<"), scaling ("<<scalesw<<"), " 
  1519. LS "smoothing ("<<smoothsw<<"), and transformation ("<<xformsw<<"). [prefix.#]" 
  1520. << ER
  1521.  
  1522. FP << leftsw << FS AP "[+int index]|'p'" AS
  1523. DP "Specifies the left cutoff index. Peak positions are reported relative" 
  1524. LS "to the left cutoff.  If 'p' is passed as the argument, the left cutoff"
  1525. LS "is set to the primer position."
  1526. << ER
  1527.  
  1528. FP << outprfxsw << FS AP "[string prefix=<input filename>]" AS
  1529. DP "Specifies file prefixes for output filenames.  If the prefix is a file,"
  1530. LS "suffixes will be added as appropriate.  If the prefix is a directory (that"
  1531. LS "is, ends in a '/'), files will be redirected to that directory and"
  1532. LS "the input filename will be used as the filename prefix.  The hash symbol"
  1533. LS "('#') may be used a placeholder for the base identifier; if it is omitted,"
  1534. LS ".# will be appended to the prefix specified with this flag."
  1535. << ER
  1536.  
  1537. FP << pmcsw << FS AP "[+double PMC=" << PMC_DEFAULT << "]" AS
  1538. DP "Specifies the minimum peak height as a product of the mean and the" 
  1539. LS "Peak Mean Coefficient (PMC).  For instance, if the mean trace value" 
  1540. LS "is 20, then -a 1.5 will only pick peaks above 24 =(1.2 * 20)." 
  1541. << ER
  1542.  
  1543. FP << prunesw << FS AP "[+int separation=" << SEPARATION_DEFAULT << "]" AS
  1544. DP "Prunes peaks by doing a pairwise comparison of adjacent peaks and"
  1545. LS "discarding the less-probably peak of each pair which is separated by less"
  1546. LS "than separation.  If os is specified, a list of pairwise comparisons which"
  1547. LS "resulted in the removal of a peak is output. [prefix.pruned]"
  1548. << ER
  1549.  
  1550. FP << pssw << FS AP
  1551. DP "Write the predicted sequence (the sequence chosen by CPeakList)" 
  1552. << ER
  1553.  
  1554. FP << ptsw << FS AP
  1555. DP "Generates traces which represent peaks and their widths & heights." 
  1556. LS "Peak traces are especially informative when overlaid with the" 
  1557. LS "individual ("<<indsw<<") trace.  [prefix.#.pt]" 
  1558. << ER
  1559.  
  1560. FP << quietsw << FS AP
  1561. << "Quiet mode.  Only errors will be displayed." 
  1562. << ER
  1563.  
  1564. FP << rightsw << FS AP "[int index]" AS
  1565. DP "Specifies the right cutoff index.  If index is <0, the right cutoff" 
  1566. LS "will be set to +index from the end of the trace.  See "<<leftsw<<"." 
  1567. << ER
  1568.  
  1569. FP << rtsw << FS AP
  1570. DP "Output the residual trace computed by subtracting the expected fluorescence"
  1571. LS "trace from the observed data [prefix.#.rt]" 
  1572. << ER
  1573.  
  1574. FP << scalesw << FS AP "[string filename]" AS
  1575. DP "Scale traces using scales from filename (in ACGT order), or default scales"
  1576. LS "if not specified." 
  1577. << ER
  1578.  
  1579. FP << smoothsw << FS AP "[short iterations="<<SMOOTH_ITER_DEFAULT<<"]" AS
  1580. DP "Smooths by using a weighted average of the 3 points about a particular"
  1581. LS "index, with the special case of the endpoints handled by throwing out the"
  1582. LS "third point and normalizing the weighting coefficients.  The process is"
  1583. LS "repeated iterations times."
  1584. << ER
  1585.  
  1586. FP << tfssw << FS AP
  1587. DP "Output a summary of tracefile statistics and peaks [prefix.tfs]" 
  1588. << ER
  1589.  
  1590. FP << tssw << FS AP
  1591. DP "Output a summary of individual trace statistics and peaks [prefix.#.ts]" 
  1592. << ER
  1593.  
  1594. FP << versionsw << FS AP
  1595. DP "Print version info" 
  1596. << ER
  1597.  
  1598. FP << vsw << FS AP
  1599. DP "(extremely) Verbose mode"
  1600. << ER
  1601.  
  1602. FP << xformsw << FS AP "[string filename]" AS
  1603. DP "Applies a 4x4 transformation/orthogonalization matrix to the 4 traces,"
  1604. LS "producing a new set of traces which replaces the existing set."
  1605. LS "Transformation matrics are expected to be 4x4 matrix of the form:"
  1606. LS "(file consists of mTS values only)"        
  1607. LS INDENT "M =         [,A]    [,C]    [,G]    [,T]"
  1608. LS INDENT "    [A,]    mAA     mAC     mAG     mAT"
  1609. LS INDENT "    [C,]    mCA     mCC     mCG     mCT"
  1610. LS INDENT "    [G,]    mGA     mGC     mGG     mGT"
  1611. LS INDENT "    [T,]    mTA     mTC     mTG     mTT"
  1612. LS "The general equation for the resulting traces is:"
  1613. LS "R = M O  <==>  R(T,i) =      sum     [ mTS x O(S,i) ]"
  1614. LS "                         S in {ACGT}"
  1615. LS "where R is the resulting vector of 4 traces"
  1616. LS "      O is the original vector of 4 traces"
  1617. LS "      M is the 4x4 ({ACGT}x{ACGT}) matrix whose elements m(i,j) are"
  1618. LS "        the cross-term contributions of channel j to channel i"
  1619. LS "      S & T are trace identifiers (Source & Target) in {A,C,G,T}"
  1620. LS "      i loops over the indices of the trace"
  1621. << ER
  1622.  
  1623. FP << zerosw << FS AP "[+double threshold=" << Z_THRESH_DEFAULT << "]" AS
  1624. DP "Specifies the epsilon about zero in which derivatives are considered to"
  1625. LS "be exactly zero, and thus the crest (trough) of a local maxima (minima)."
  1626. << ER;
  1627.  
  1628. cout <<
  1629. LS "Argument specifications are given in square brackets ('[' & ']') with their"
  1630. LS "associated flags where applicable.  Usually exactly one argument is"
  1631. LS "expected; 0+ means that zero or more arguments may follow; 0,1 indicates"
  1632. LS "0 or 1 arguments are expected (ie. it's optional).  If a flag has an"
  1633. LS "optional argument which is omitted by the user, the following flag must"
  1634. LS "begin with a hyphen ('-'); otherwise, the hyphen is optional.  Default"
  1635. LS "values are designated in the argument specification by an equal sign ('=')."
  1636. LS "If any operation (ie. a flag that results in an output file) is specified,"
  1637. LS "then none of the default operations are performed unless specifically"
  1638. LS "requested.  If a flag causes an output file to be written, the default"
  1639. LS "filename for that file is indicated in square brackets.  Generally,"
  1640. LS "filnames are generated by appending a suffix to the input filename and are"
  1641. LS "written to the same directory as the input file; the user may specify an"
  1642. LS "alternate prefix or destination directory with the -o flag.  The hash"
  1643. LS "symbol ('#') is a placeholder for the base letter when the output file"
  1644. LS "writes base-specific information; # may be used on the command line.  If"
  1645. LS "filename ends with '.s1' the suffix is removed.  Peaks will be picked"
  1646. LS "only if necessary for the requested operations."
  1647. LS "WARNING: If file with an output filename already exists, it will be"
  1648. LS "         overwritten.";
  1649.     }
  1650.